home *** CD-ROM | disk | FTP | other *** search
/ GFX Sensations 1 / Graphic Sensations - Volume 1.iso / tools / amiga / 3d_tools / irit40s.lha / Irit / geom_lib / convex.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-30  |  16.8 KB  |  508 lines
  1. /******************************************************************************
  2. * Convex.c - test convexity and converts polygons to convex ones.          *
  3. *******************************************************************************
  4. * Written by Gershon Elber, March 1990.                          *
  5. ******************************************************************************/
  6.  
  7. /* #define DEBUG2           Defines some printing/debugging routines. */
  8.  
  9. #include <stdio.h>
  10. #include <ctype.h>
  11. #include <math.h>
  12. #include <string.h>
  13. #include "allocate.h"
  14. #include "convex.h"
  15. #include "poly_cln.h"
  16. #include "geomat3d.h"
  17. #include "intrnrml.h"
  18. #include "priorque.h"
  19.  
  20. /* Used to hold edges (V | V -> Pnext) that intersect with level y = Ylevel. */
  21. typedef struct InterYVrtxList {
  22.     ByteType InterYType;
  23.     IPVertexStruct *V;
  24.     struct InterYVrtxList *Pnext;
  25. } InterYVrtxList;
  26.  
  27. #define CONVEX_EPSILON  1e-4       /* Colinearity of two normalized vectors. */
  28.  
  29. #define    INTER_Y_NONE    0               /* Y level intersection type. */
  30. #define    INTER_Y_START    1
  31. #define    INTER_Y_MIDDLE    2
  32.  
  33. #define LOOP_ABOVE_Y    0      /* Type of open loops extracted from polygon. */
  34. #define LOOP_BELOW_Y    1
  35.  
  36. /* Used to sort and combine the polygons above Ylevel together if possible.  */
  37. typedef struct SortPolysInX {
  38.     IPVertexStruct *VMinX, *VMaxX;
  39.     int Reverse;      /* If TRUE than VMinX vertex is AFTER VMaxX vertex. */
  40. } SortPolysInX;
  41.  
  42. /* The following are temporary flags used to mark vertices that were visited */
  43. /* by the loop tracing, at list once. As each vertex may be visited two at   */
  44. /* the most (as starting & as end point of open loop), this is enough.         */
  45. /* INTER_TAG is used to mark vertices that created brand new with intersected*/
  46. /* with the line y = Ylevel. Those are used to detect INTERNAL edges - if    */
  47. /* at list one end of it is INTER_TAG, that edge is INTERNAL.             */
  48. #define INTER_TAG   0x40
  49. #define VISITED_TAG 0x80
  50.  
  51. #define    IS_INTER_VRTX(Vrtx)    ((Vrtx)->Tags & INTER_TAG)
  52. #define    SET_INTER_VRTX(Vrtx)    ((Vrtx)->Tags |= INTER_TAG)
  53. #define    RST_INTER_VRTX(Vrtx)    ((Vrtx)->Tags &= ~INTER_TAG)
  54.  
  55. #define    IS_VISITED_VRTX(Vrtx)    ((Vrtx)->Tags & VISITED_TAG)
  56. #define    SET_VISITED_VRTX(Vrtx)    ((Vrtx)->Tags |= VISITED_TAG)
  57. #define    RST_VISITED_VRTX(Vrtx)    ((Vrtx)->Tags &= ~VISITED_TAG)
  58.  
  59. static int SplitPolyIntoTwo(IPPolygonStruct *Pl, IPVertexStruct *V,
  60.                 IPPolygonStruct **Pl1, IPPolygonStruct **Pl2);
  61. static IPVertexStruct *FindRayPolyInter(IPPolygonStruct *Pl,
  62.         IPVertexStruct *VRay, PointType RayDir, PointType PInter);
  63. static void TestConvexityDir(IPPolygonStruct *Pl);
  64.  
  65. #ifdef DEBUG2
  66. static void PrintVrtxList(IPVertexStruct *V);
  67. static void PrintPoly(IPPolygonStruct *P);
  68. #endif /* DEBUG2 */
  69.  
  70. /*****************************************************************************
  71. *   Routine to prepare transformation martix to rotate such that Dir is         *
  72. * parallel to the Z axes. Used by the convex decomposition to rotate the     *
  73. * polygons to be XY plane parallel.                         *
  74. *    Algorithm: form a 4 by 4 matrix from Dir as follows:             *
  75. *                |  Tx  Ty  Tz  0 |   A transformation which takes the coord *
  76. *                |  Bx  By  Bz  0 |  system into T, N & B as required.         *
  77. * [X  Y  Z  1] * |  Nx  Ny  Nz  0 |                         *
  78. *                |  0   0   0   1 |                         *
  79. * N is exactly Dir, but we got freedom on T & B - T & B must be on         *
  80. * a plane perpendicular to N and perpendicular between them but thats all!   *
  81. * T is therefore selected using this (heuristic ?) algorithm:             *
  82. * Let P be the axis of which the absolute N coefficient is the smallest.     *
  83. * Let B be (N cross P) and T be (B cross N).                     *
  84. *****************************************************************************/
  85. void GenRotateMatrix(MatrixType Mat, VectorType Dir)
  86. {
  87.     int i, j;
  88.     RealType R;
  89.     VectorType DirN, T, B, P;
  90.  
  91.     PT_COPY(DirN, Dir);
  92.     PT_NORMALIZE(DirN);
  93.     PT_CLEAR(P);
  94.     for (i = 1, j = 0, R = ABS(DirN[0]); i < 3; i++)
  95.     if (R > ABS(DirN[i])) {
  96.         R = DirN[i];
  97.         j = i;
  98.     }
  99.     P[j] = 1.0;/* Now P is set to the axis with the biggest angle from DirN. */
  100.  
  101.     GMVecCrossProd(B, DirN, P);                  /* calc the bi-normal. */
  102.     PT_NORMALIZE(B);
  103.     GMVecCrossProd(T, B, DirN);                /* calc the tangent. */
  104.  
  105.     MatGenUnitMat(Mat);
  106.     for (i = 0; i < 3; i++) {
  107.     Mat[i][0] = T[i];
  108.     Mat[i][1] = B[i];
  109.     Mat[i][2] = DirN[i];
  110.     }
  111. }
  112.  
  113. /*****************************************************************************
  114. *   Routine to test all polygons in a given object for convexity, and split  *
  115. * non convex ones, non destructively - the original object is not modified.  *
  116. *****************************************************************************/
  117. IPObjectStruct *ConvexPolyObjectN(IPObjectStruct *PObj)
  118. {
  119.     IPObjectStruct
  120.     *PObjCopy= CopyObject(NULL, PObj, FALSE);
  121.  
  122.     ConvexPolyObject(PObjCopy);
  123.  
  124.  
  125.     return PObjCopy;
  126. }
  127.  
  128. /*****************************************************************************
  129. *   Routine to test all the polygons in a given object for convexity, and    *
  130. * split the non convex ones.                             *
  131. *****************************************************************************/
  132. void ConvexPolyObject(IPObjectStruct *PObj)
  133. {
  134.     IPPolygonStruct *Pl, *PlSplit, *PlTemp,
  135.     *PlPrev = NULL;
  136.  
  137.     if (!IP_IS_POLY_OBJ(PObj))
  138.     IritFatalError("ConvexPolyObject: parameter given is not polygonal object");
  139.  
  140.     if (IP_IS_POLYLINE_OBJ(PObj)) {
  141.     return;
  142.     }
  143.  
  144.     Pl = PObj -> U.Pl;
  145.     while (Pl != NULL) {
  146.     if (!ConvexPolygon(Pl)) {
  147.         PlSplit = SplitNonConvexPoly(Pl);
  148.         CleanUpPolygonList(&PlSplit);       /* Zero length edges etc. */
  149.         if (PlSplit != NULL) {     /* Something is wrong here, ignore. */
  150.         if (Pl == PObj -> U.Pl)
  151.             PObj -> U.Pl = PlSplit;               /* First. */
  152.         else
  153.             PlPrev -> Pnext = PlSplit;
  154.  
  155.         PlTemp = PlSplit;
  156.         while (PlTemp -> Pnext != NULL)
  157.             PlTemp = PlTemp -> Pnext;
  158.         PlTemp -> Pnext = Pl -> Pnext;
  159.         PlPrev = PlTemp;
  160.         IPFreePolygon(Pl);            /* Free old polygon. */
  161.         Pl = PlPrev -> Pnext;
  162.         }
  163.         else {
  164.         if (Pl == PObj -> U.Pl)
  165.             PObj -> U.Pl = Pl -> Pnext;
  166.         PlPrev = Pl -> Pnext;
  167.         IPFreePolygon(Pl);            /* Free old polygon. */
  168.         Pl = PlPrev;
  169.         }
  170.     }
  171.     else {
  172.         PlPrev = Pl;
  173.         Pl = Pl -> Pnext;
  174.     }
  175.     }
  176. }
  177.  
  178. /*****************************************************************************
  179. *   Routine to test if the given polygon is convex or not.             *
  180. * Algorithm: The polygon is convex iff the normals generated from cross      *
  181. * products of two consecutive edges points to the same direction. he same    *
  182. * direction is tested by dot product with the polygon plane normal which     *
  183. * should produce consistent sign for all normals, in order the polygon to    *
  184. * be convex.                                     *
  185. *   The routine returns TRUE iff the polygon is convex. In addition the      *
  186. * polygon CONVEX tag (see IPPolygonStruct) is also updated.             *
  187. *   Note that if the polygon IS marked as convex, nothing is tested!         *
  188. *   Also convex polygons are tested so that the vertices are ordered such    *
  189. * that cross product of two adjacent edges will be in the normal direction.  *
  190. *****************************************************************************/
  191. int ConvexPolygon(IPPolygonStruct *Pl)
  192. {
  193.     int FirstTime = TRUE;
  194.     RealType Size,
  195.     NormalSign = 0.0;
  196.     VectorType V1, V2, PolyNormal, Normal;
  197.     IPVertexStruct *VNext,
  198.     *V = Pl -> PVertex;
  199.  
  200.     if (IP_IS_CONVEX_POLY(Pl))
  201.     return TRUE;                   /* Nothing to do around here. */
  202.  
  203.     /* Copy only A, B, C from Ax+By+Cz+D = 0: */
  204.     PT_COPY(PolyNormal, Pl -> Plane);
  205.  
  206.     do {
  207.     VNext = V -> Pnext;
  208.  
  209.     PT_SUB(V1, VNext -> Coord, V -> Coord);
  210.     if ((Size = PT_LENGTH(V1)) > EPSILON) {
  211.         Size = 1.0 / Size;
  212.         PT_SCALE(V1, Size);
  213.     }
  214.     PT_SUB(V2, VNext -> Pnext -> Coord, VNext -> Coord);
  215.     if ((Size = PT_LENGTH(V2)) > EPSILON) {
  216.         Size = 1.0 / Size;
  217.         PT_SCALE(V2, Size);
  218.     }
  219.     GMVecCrossProd(Normal, V1, V2);
  220.  
  221.     if (PT_LENGTH(Normal) < CONVEX_EPSILON) {
  222.         V = VNext;
  223.         continue;                    /* Skip colinear points. */
  224.     }
  225.  
  226.     if (FirstTime) {
  227.         FirstTime = FALSE;
  228.         NormalSign = DOT_PROD(Normal, PolyNormal);
  229.     }
  230.     else if (NormalSign * DOT_PROD(Normal, PolyNormal) < 0.0) {
  231.         IP_RST_CONVEX_POLY(Pl);
  232.         return FALSE;          /* Different signs --> not convex. */
  233.     }
  234.  
  235.     V = VNext;
  236.     }
  237.     while (V != Pl -> PVertex);
  238.  
  239.     IP_SET_CONVEX_POLY(Pl);
  240.  
  241.     if (NormalSign < 0.0)
  242.         IritPrsrReverseVrtxList(Pl);
  243.  
  244.     return TRUE;    /* All signs are the same --> the polygon is convex. */
  245. }
  246.  
  247. /*****************************************************************************
  248. *   Routine to split non convex polygon to list of convex ones.             *
  249. * This algorithm is much simpler than the one used before:             *
  250. * 1. Remove a polygon from GlblList. If non exists stop.             *
  251. * 2. Search for non convex corner. If not found stop - polygon is convex.    *
  252. *    Otherwise let the non convex polygon found be V(i).             *
  253. * 3. Fire a ray from V(i) in the opposite direction to V(i-1). Find the      *
  254. *    closest intersection of the ray with polygon boundary P.             *
  255. * 4. Split the polygon into two at V(i)-P edge and push the two new polygons *
  256. *    on the GlblList.                                 *
  257. * 5. Goto 1.                                     *
  258. *****************************************************************************/
  259. IPPolygonStruct *SplitNonConvexPoly(IPPolygonStruct *Pl)
  260. {
  261.     int IsConvex;
  262.     RealType Size;
  263.     IPPolygonStruct *GlblList, *Pl1, *Pl2,
  264.     *GlblSplitPl = NULL;
  265.     VectorType V1, V2, PolyNormal, Normal;
  266.     IPVertexStruct *V, *VNext;
  267.  
  268.     TestConvexityDir(Pl);
  269.  
  270.     GlblList = IPAllocPolygon(0, 0, CopyVertexList(Pl -> PVertex), NULL);
  271.     PLANE_COPY(GlblList -> Plane, Pl -> Plane);
  272.  
  273.     /* Copy only A, B, C from Ax+By+Cz+D = 0 plane equation: */
  274.     PT_COPY(PolyNormal, Pl -> Plane);
  275.  
  276.     while (GlblList != NULL) {
  277.     Pl = GlblList;
  278.     GlblList = GlblList -> Pnext;
  279.     Pl -> Pnext = NULL;
  280.  
  281.     IsConvex = TRUE;
  282.     V = Pl -> PVertex;
  283.     do {
  284.         VNext = V -> Pnext;
  285.  
  286.         PT_SUB(V1, VNext -> Coord, V -> Coord);
  287.          if ((Size = PT_LENGTH(V1)) > EPSILON) {
  288.         Size = 1.0 / Size;
  289.         PT_SCALE(V1, Size);
  290.         }
  291.         PT_SUB(V2, VNext -> Pnext -> Coord, VNext -> Coord);
  292.          if ((Size = PT_LENGTH(V2)) > EPSILON) {
  293.         Size = 1.0 / Size;
  294.         PT_SCALE(V2, Size);
  295.         }
  296.         GMVecCrossProd(Normal, V1, V2);
  297.         if (PT_LENGTH(Normal) < CONVEX_EPSILON) {
  298.         V = VNext;
  299.         continue;                /* Skip colinear points. */
  300.         }
  301.  
  302.         if (DOT_PROD(Normal, PolyNormal) < 0.0 &&
  303.         SplitPolyIntoTwo(Pl, V, &Pl1, &Pl2)) {
  304.         Pl -> PVertex = NULL; /* Dont free vertices - used in Pl1/2. */
  305.         IPFreePolygon(Pl);
  306.  
  307.         Pl1 -> Pnext = GlblList;       /* Push polygons on GlblList. */
  308.         GlblList = Pl1;
  309.         Pl2 -> Pnext = GlblList;
  310.         GlblList = Pl2;
  311.  
  312.         IsConvex = FALSE;
  313.         break;
  314.         }
  315.  
  316.         V = VNext;
  317.     }
  318.     while (V != Pl -> PVertex);
  319.  
  320.     if (IsConvex) {
  321.         IP_SET_CONVEX_POLY(Pl);
  322.         Pl -> Pnext = GlblSplitPl;
  323.         GlblSplitPl = Pl;
  324.     }
  325.     }
  326.  
  327.     return GlblSplitPl;
  328. }
  329.  
  330. /*****************************************************************************
  331. *   Split polygon at the vertex specified by V -> Pnext, given V, into two,  *
  332. * by firing a ray from V -> Pnext in the opposite direction to V and finding *
  333. * closest intersection with the polygon P. (V -> Pnext, P) is the edge to    *
  334. * split the polygon at.                                 *
  335. *****************************************************************************/
  336. static int SplitPolyIntoTwo(IPPolygonStruct *Pl, IPVertexStruct *V,
  337.                 IPPolygonStruct **Pl1, IPPolygonStruct **Pl2)
  338. {
  339.     PointType Vl, PInter;
  340.     IPVertexStruct *VInter, *VNew1, *VNew2;
  341.  
  342.     PT_SUB(Vl, V -> Pnext -> Coord, V -> Coord);
  343.     VInter = FindRayPolyInter(Pl, V -> Pnext, Vl, PInter);
  344.     V = V -> Pnext;
  345.  
  346.     if (VInter == NULL ||
  347.     VInter == V ||
  348.     VInter -> Pnext == V)
  349.     return FALSE;
  350.  
  351.     /* Make the two polygon vertices lists. */
  352.     VNew1 = IPAllocVertex(V -> Count, V -> Tags, NULL, V -> Pnext);
  353.     PT_COPY(VNew1 -> Coord, V -> Coord);
  354.     PT_COPY(VNew1 -> Normal, V -> Normal);
  355.     IP_SET_INTERNAL_VRTX(V);
  356.     if (PT_APX_EQ(VInter -> Coord, PInter)) {
  357.     /* Intersection points is close to VInter point. */
  358.     VNew2 = IPAllocVertex(VInter -> Count, VInter -> Tags, NULL,
  359.                             VInter -> Pnext);
  360.     PT_COPY(VNew2 -> Coord, VInter -> Coord);
  361.     PT_COPY(VNew2 -> Normal, VInter -> Normal);
  362.     VInter -> Pnext = VNew1;
  363.     IP_SET_INTERNAL_VRTX(VInter);
  364.     V -> Pnext = VNew2;
  365.     }
  366.     else if (PT_APX_EQ(VInter -> Pnext -> Coord, PInter)) {
  367.     /* Intersection points is close to VInter -> Pnext point. */
  368.     VNew2 = IPAllocVertex(VInter -> Pnext -> Count,
  369.                   VInter -> Pnext -> Tags,
  370.                   NULL, VInter -> Pnext -> Pnext);
  371.     PT_COPY(VNew2 -> Coord, VInter -> Pnext -> Coord);
  372.     PT_COPY(VNew2 -> Normal, VInter -> Pnext -> Normal);
  373.     VInter -> Pnext -> Pnext = VNew1;
  374.     IP_SET_INTERNAL_VRTX(VInter -> Pnext);
  375.     V -> Pnext = VNew2;
  376.     }
  377.     else {
  378.     /* PInter is in the middle of (VInter, VInter -> Pnext) edge: */
  379.     VNew2 = IPAllocVertex(VInter -> Count, VInter -> Tags, NULL,
  380.                             VInter -> Pnext);
  381.     PT_COPY(VNew2 -> Coord, PInter);
  382.     VInter -> Pnext = IPAllocVertex(0, 0, NULL, VNew1);
  383.     PT_COPY(VInter -> Pnext -> Coord, PInter);
  384.     InterpNrmlBetweenTwo(VNew2, VInter, VNew2 -> Pnext);
  385.     PT_COPY(VInter -> Pnext -> Normal, VNew2 -> Normal);
  386.     IP_SET_INTERNAL_VRTX(VInter -> Pnext);
  387.     V -> Pnext = VNew2;
  388.     }
  389.  
  390.     *Pl1 = IPAllocPolygon(0, 0, VNew1, NULL);
  391.     PLANE_COPY((*Pl1) -> Plane, Pl -> Plane);
  392.     *Pl2 = IPAllocPolygon(0, 0, VNew2, NULL);
  393.     PLANE_COPY((*Pl2) -> Plane, Pl -> Plane);
  394.  
  395.     return TRUE;
  396. }
  397.  
  398. /*****************************************************************************
  399. *   Find where a ray first intersect a given polygon. The ray starts at one  *
  400. * of the polygon vertices and so distance less than EPSILON is ignored.      *
  401. *   Returns the vertex V in which (V, V -> Pnext) has the closest         *
  402. * intersection with the ray PRay, DRay at Inter.                 *
  403. *****************************************************************************/
  404. static IPVertexStruct *FindRayPolyInter(IPPolygonStruct *Pl,
  405.             IPVertexStruct *VRay, PointType RayDir, PointType PInter)
  406. {
  407.     RealType t1, t2,
  408.     MinT = INFINITY;
  409.     PointType Vl, Ptemp1, Ptemp2, PRay;
  410.     IPVertexStruct *VNext,
  411.     *V = Pl -> PVertex,
  412.     *VInter = NULL;
  413.  
  414.     PT_COPY(PRay, VRay -> Coord);
  415.  
  416.     do {
  417.     VNext = V -> Pnext;
  418.     if (V != VRay && VNext != VRay) {
  419.         PT_SUB(Vl, VNext -> Coord, V -> Coord);
  420.         if (CGDistPointLine(PRay, V -> Coord, Vl) > EPSILON) {
  421.         /* Only if the point the ray is shoot from is not on line: */
  422.         CG2PointsFromLineLine(PRay, RayDir, V -> Coord, Vl,
  423.                       Ptemp1, &t1, Ptemp2, &t2);
  424.         if (CGDistPointPoint(Ptemp1, Ptemp2) < EPSILON * 10.0 &&
  425.             t1 > EPSILON && t1 < MinT &&
  426.             (t2 <= 1.0 || APX_EQ(t2, 1.0)) &&
  427.             (t2 >= 0.0 || APX_EQ(t2, 0.0))) {
  428.             PT_COPY(PInter, Ptemp2);
  429.             VInter = V;
  430.             MinT = t1;
  431.         }
  432.         }
  433.     }
  434.     V = VNext;
  435.     }
  436.     while (V != Pl -> PVertex && V -> Pnext != NULL);
  437.  
  438.     return VInter;
  439. }
  440.  
  441. /*****************************************************************************
  442. *   Test convexity direction - a cross product of two edges of a convex      *
  443. * corner of the polygon should point in the normal direction. if this is not *
  444. * the case - the polygon vertices are reveresed.                 *
  445. *****************************************************************************/
  446. static void TestConvexityDir(IPPolygonStruct *Pl)
  447. {
  448.     int Coord = 0;
  449.     VectorType V1, V2, Normal;
  450.     IPVertexStruct *V, *VExtrem;
  451.  
  452.     /* Find the minimum component direction of the normal and used that axes */
  453.     /* to find an extremum point on the polygon - that extrmum point must be */
  454.     /* a convex corner so we can find the normal direction of convex corner. */
  455.     if (ABS(Pl -> Plane[1]) < ABS(Pl -> Plane[Coord]))
  456.     Coord = 1;
  457.     if (ABS(Pl -> Plane[2]) < ABS(Pl -> Plane[Coord]))
  458.     Coord = 2;
  459.     V = VExtrem = Pl -> PVertex;
  460.     do {
  461.     if (V -> Coord[Coord] > VExtrem -> Coord[Coord])
  462.         VExtrem = V;
  463.     V = V -> Pnext;
  464.     }
  465.     while (V != Pl -> PVertex && V != NULL);
  466.  
  467.     /* Make sure next vertex is not at the extremum value: */
  468.     while (APX_EQ(VExtrem -> Coord[Coord], VExtrem -> Pnext -> Coord[Coord]))
  469.     VExtrem = VExtrem -> Pnext;
  470.  
  471.     /* O.K. V form a convex corner - evaluate its two edges cross product:   */
  472.     for (V = Pl -> PVertex; V -> Pnext != VExtrem; V = V -> Pnext); /* Prev. */
  473.     PT_SUB(V1, VExtrem -> Coord, V -> Coord);
  474.     PT_SUB(V2, VExtrem -> Pnext -> Coord, VExtrem -> Coord);
  475.     GMVecCrossProd(Normal, V1, V2);
  476.  
  477.     if (DOT_PROD(Normal, Pl -> Plane) < 0.0)
  478.     IritPrsrReverseVrtxList(Pl);
  479. }
  480.  
  481. #ifdef DEBUG2
  482.  
  483. /*****************************************************************************
  484. *   Print the content of the given vertex list, to standard output.         *
  485. *****************************************************************************/
  486. static void PrintPoly(IPPolygonStruct *P)
  487. {
  488.     IPVertexStruct
  489.     *VFirst = P -> PVertex,
  490.     *V = VFirst;
  491.  
  492.     if (V == NULL)
  493.     return;
  494.  
  495.     printf("VERTEX LIST:\n");
  496.     do {
  497.     printf("%12lg %12lg %12lg, Tags = %02x\n",
  498.         V -> Coord[0], V -> Coord[1], V -> Coord[2], V -> Tags);
  499.     V = V -> Pnext;
  500.     }
  501.     while (V != NULL && V != VFirst);
  502.  
  503.     if (V == NULL)
  504.     printf("Loop is not closed (NULL terminated)\n");
  505. }
  506.  
  507. #endif /* DEBUG2 */
  508.